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I. INTRODUCTION 


iweb) is an integrable function of the real variable 


t, the function f(s) defined by the integral 


со 


f(s) -| e 8* r(t)gt (1) 


O 


is called the Laplace transformation, or the Laplace trans- 


form, of F(t) and is frequently indicated by the notation 
its) = [F(t)] (2) 


The Laplace transformation is a linear integral trans- 
formation which is widely used in applied mathematics and 
technology. A typical application is as follows. An un- 
known function F(t) satisfies a certain differential equa- 
tion and specified end conditions. Employing theorems 
applicable to the Laplace transform, a function f(s) is 
determined as the solution of an algebraic equation which 
corresponds, in an appropriate fashion, to the differential 
equation and end conditions of F(t). Then the problem is 
reduced to the following task: having established f(s), 


it is necessary to determine the function F(t) for which 


aes) 23 [F(t)]. 





This problem is called finding the inverse (Laplace) 
Mans torm of the function f(s). This is customarily 


written as 
F(t) = *{£(s)] (3) 


It may be shown that the following inversion integral 


F(t) = == J её f(s)ds (4) 


C 


accomplishes the desired inversion. Jn this expression 


s is the complex variable 
S = xXx tiy СЭ) 


апа the integration is along a suitable path in the complex 
s-plane. 

п tact there is a function F(t) such that the given 
f(s) = 3 [F(t)], then the function determined by equation 
(4) is, indeed, the function F(t). The result is unique 


in the following sense. If 


f(s) -| CE uch (6) 
@ 
апа 
F,(t) = == £(s)ds (7) 





then F. (t) and Fo(t) differ only on a set of Lebesgue 
measure zero. 

In this thesis numerical methods will be utilized to 
perform the integration of equation (4) and functions for 
which the inversion is successful will correspond to 
functions F(t) which are at least piecewise continuous, 
having only isolated discontinuities, if any are present 
at all. If F(t) is discontinuous for t = a, then the 


recovered function satisfies the condition 


F(a) = 5 Lim[F(ate) + F(a-e)] (8) 


ЕО 


The path of integration, indicated by C in equation 
(4) is frequently described as an infinitely long, vertical 
straight line x = constant, to the right of any singulari- 
ties of f(s). (A theorem of complex variables shows that 
f(s) is analytic in its region of convergence, a right-hand 
plane, except for isolated singularities.) This straight 
vertical contour is frequently called the Bromwich contour. 

Usually inversion of a Laplace transform, that is find- 
ime F(t) eee | es) is accomplished by use of a table of 
transform pairs. One of the most complete of these is that 
of Roberts and Kaufman [14]. Other such tables are listed 
@sereterences [5], [6], [7], and [13]. 

several theorems concerning the Laplace transform extend 
the usefulness of such transformation pair tables. For 


example 


10 





Chie (s) + 2,(8)] = [21 (s)] +O te, (8)! 
(9) 
However, it is not an infrequent occurrence that one is 
called upon to find F(t) in a case where the tables and 
the theorems are of no direct help. There has developed 
a considerable literature concerning this problem. 
One obvious method is to attempt to expand the given 
ШЕ аз a series of functions f. Cs) for which the inverses 
are known. Then using the linear property indicated in 


equation (9), the result may be expressed as 
s -l(e(s)] -A lísa f. (s)] = sa £ l[# (s)] (10) 


The integration indicated in equation (4) can also be 
performed by various analytical and approximate methods. 
In particular, one can perform the integration along the 
Bromwich contour by a numerical procedure. Much of the 
literature is devoted to such methods. 

In general these methods amount to the following. А 
finite number of points along the Bromwich contour are 
determined according to some law and the integrand is 
evaluated at such points. There are associated weight fac- 
tors and the integral is approximated as a sum of products 
of integrand values times weight factors. One of the most 
widely used of such methods is described by Salzer [15], 
[16], and [17]. Other procedures of the type are treated 


ИШ ОШО 120], [2], [8], [9], [10], and [21]. 


JI 





It is interesting to note that the same classical poly- 
nomials which are useful in the expansions indicated by 
equation (10) are also encountered in locating points for 
evaluation using a numerical integration scheme so that 
the final algorithm may be the same even though the original 
motivation was quite different. 

Other methods which do not fall clearly into the 
classes of methods discussed above have been presented by 
Widder [22], and Bellman, Kalaba, and Lockett [3]. 

These approximate methods have not always been success- 
ful. Hiep [12] employed Salzer's method for numerical in- 
version on a problem of heat transfer in porous media and 
und it difficult to obtain results which did not exhibit 
physically incorrect behavior, e.g., in the vicinity of a 
point of engineering interest, the temperature of the 
medium fell below sink temperature and rose above source 
temperature. In similar work on a conjugated heat transfer 
problem Zargary [23] encountered the same difficulty, and, 
after devoting considerable time and attention to the prob- 
lem, was forced to abandon the use of numerical inversion. 

To overcome such difficulties in problems of these 
kinds, we have reverted to direct numerical integration of 
equation (4) along appropriate contours in the s-plane. 

Our procedures, which appear thus far to be efficient and 


reliable, are described in what follows. 


12 





IIS DI STOSTED CONTOURS 


We have found that great advantage accrues from distort- 
ing the Bromwich contour in ways which will subsequently 
be described. 

According to theorems concerning analytic functions of 
a complex variable, the Bromwich contour can be distorted 
arbitrarily as long as the process of distortion does not 
cause the contour C to cross over a Singularity of f(s) and 
as long as the two ends of C extend to infinity in a manner 
which preserves the convergence of the integral in equation 
(4). 

Hides. with the supposition that f(s) is analytic, with 
the exception of isolated singularities, and is real valued 


for real s, let g(s) represent the integrand 


g(s) = et f(s) = u(s) + iv(s) Ob 


where u(s) and v(s) are the real and imaginary harmonic 
components of g(s). Furthermore, assume that the contour C 
consists of the two symmetrical parts с” апа c as shown 

in Figure 1. Consideration of the contributions PORT [f(s)] 


* жж 
over arcs ds, On © and ds, on C yields the following: 


dx, dx, 
dy, = dy; (12) 
т 
Vo — 
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FIGURE 1. Distorted Bromwich Contour. Contour is. 
distorted into C* and its reflection C . 


caus, 
I — a 
511181981 + gods,] = —[vydx y + uydyj] se) 


and, therefore, 


ЕСЕ) | Ge eerily | (14) 


x 


e 


u 
ale 


САША ао (14) is, of course, further simplified if on 
is selected to be a path along which v=0. Additional 
theoretical reasons exist for such a selection. Such a path 
will subsequently be referred to as a path of steepest 
descent. Numerical integration along such a path, ora 
close facsimile thereof, was done by Esch [11] in 1957 and 


Carrier, Krook, and Pearson [4] recommended further 


14 





exploration of this procedure. However, no further such in- 
vestigation seems to have been done until the present thesis. 
In Chapter V, a description of an algorithm for performing 
numerical integration along a path oo for which v=0 over 
a substantial portion of its range, will be presented. 
Almost by accidental discovery, the present work has 
revealed that it is computationally more efficient simply 
to make an arbitrary choice for the contour с” which satis- 
fies the following criteria: 


* 
(1) C is a smooth continuous curve of the form 


eo Хр + iy tp) (T9) 


where p is a real parameter O < p < °. 
(2) с“ lies to the right of all Singularities of f(s). 
(ES) в. does not extend sufficiently far to the right 
to cause computational difficulties due to the 
factor . 
(4) с” approaches infinity in such a way that x 
approaches -~ and y approaches œ, 
(5) Oscillation of the integrand is not excessive. 
(Otherwise there may be loss of accuracy due to 
the sampling rate of the integration algorithm 
seriously mismatching variations in the integrand 
and/or positive and negative contributions nearly 
cancelling each other.) 
Such a contour will be referred to as a Simple parameterized 


contour. Discussion of how to assure satisfying criteria 


15 





(2), (3), and (5) will be deferred until after the develop- 
Ment Of the actual algorithm. 
The integral shown in equation (14) is approximated by 


the simple trapezoidal sum 


M 
P(t) = gz 2 (ives tvs, 11% ху] 
k=] 


+ [u(s, tu Sk) Yke Ygl? (16) 


where 


S. “= X. + ly, = x€p,)+ iy(p,) (17) 


The arbitrarily chosen functions x(p) and y(p) are monotoni- 


cally increasing and Py» Por+++, P 1S an increasing 


n 
sequence. 


` 


A very simple example, but one which has been employed 


very successfully, is given by the following equations 


x = А - Bp? 


Ye) (18) 


p = 0, А О ee 


eat number A is selected so as to satisfy criteria (2) 
and (3) above. The real number B has, usually, been 
assigned a value of one. 

The sum includes only a finite number of terms. It is 
necessary to provide an appropriate criterion for terminat- 


ing the process of summation. This will be discussed later. 


16 





Mmonaere tO monitor the possibility of inaccurate 
evaluations resulting from oscillations,in addition to the 


sum F(t) given by equation (16), the sum G(t) given by 


M 
000) = 9 УХ vC ls x] 
к= 1 
+ [uts,,„)tucs )] X yasa) (19) 


is also obtained. This is the sum of the absolute values 
of the addends to the sum (16). Additionally, the number 
of times successive addends are of opposite signs is 
recorded. 

The present numerical experiments, to be described 
later, were generally so successful that the effect of 
oscillation could not be seen. However, in some cases 
where poor contours с” were deliberately chosen, it was 
found that inaccurate results were obtained if G(t) was 


many orders of magnitude (e.g., 10°) times as great as F(t). 
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III. THE COMPUTER PROGRAM AND 
ITS PRINCIPAL SUBROUTINES 

Previous discussion proposed two alternatives in terms 
of advantageous distortion of the Bromwich contour upon 
Which to perform the numerical integration of equation (16). 
The simple parameterized curve has been experimentally found 
to be the more efficient of the two methods and its algo- 
rithm will be developed at this time. The steepest descent 
procedure will be treated in Chapter V. 

The computer program and subprograms which implement 
the simple procedure developed in the preceding chapter are 
all written in the FORTRAN language, using double precision 
arithmetic. They have been tested and debugged on an IBM 
360/67. For machines having a larger mantissa, single pre- 
F n; arithmetic might prove to be satisfactory. 

SUBROUTINE VALUE is a user supplied routine which imple- 
ments computation of the g(s) given by equation (11) for the 
desired transform f(s). This subroutine calculates the u(s) 
and v(s) corresponding to a given x and y. 

Either complex or real arithmetic may be utilized in this 
program, although complex arithmetic offers a decided advan- 
tage in convenience and simplicity. However, the user 
Should be aware of the potential dangers in unintentional 
exchange of sheets of a Riemann surface when investigating 


a multi-valued function, More will be said about this later. 
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Also, although FORTRAN manuals state or imply that the 
FUNCTIONS DREAL and DIMAG are implemented in all compilers, 
this may not actually be the case. Accordingly it is prudent 
to add the following statements after SUBROUTINE VALUE. 

FUNCTION DREAL(CPLX) 
REAL*S CPLX(2) 

DREAL = CPLX(1) 
RETURN 

END 

FUNCTION DIMAG(CPLX) 
REAL*S CPLX(2) 

DIMAG = CPLX(2) 
RETURN 

END 

SUBROUTINE CURVE is also a user supplied routine which 
calculates the x and y values for each consecutive incrementa- 
meen Of the parameter p of equation (17), thus locating a 
point on the simple parameterized contour. Equation (18) 
provides an illustrative and highly effective example of 
parabolic form which has been employed almost exclusively 
during this investigative work. 


SUBROUTINE INCREP formulates the increasing sequence of 


the parameter p as given by the following equation 


Uude ( 20) 


G 
il 
jo) 


Limited analysis with other than equally spaced incrementa- 
Ono points along the contour of integration did not ex- 
hibit any enhancing capabilities and was not continued. 

The main program accomplishes the integration of equation 
(14) using a trapezoidal summation as described by equation 


EU In addition, the absolute value of the integrand is 
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summed in accordance with equation (19). Also, the number 
of times a sign changes occurs between successive evaluations 
Of the addend of equation (16) is recorded. These features 
afford the user an opportunity to monitor potential oscilla- 
tion effects upon the accuracy of the results. Clock time 
and the number of calls to SUBROUTINE VALUE required for 
performing the numerical contour integration are maintained 
as measures of the computational efficiency of the numerical 
inversion procedure. 

Termination of numerical integration occurs when each 
Of N successive evaluations of the addend of equation (19) 
has magnitude less than a specified epsilon. We have usually 
used N = 5, This is a very stringent requirement and could 
probably be relaxed resulting in some saving of computational 
time and with but minor loss of accuracy. 

In addition to the termination criterion described in 
the previous paragraph, other program inputs include the 
increment A of equation (13), appropriate numerical values 
of any constants of the function f(s), the starting position 
A of equation (18) and the value of t. 

fee program output prints the starting position of the 
integration contour, the values given by equations (16) and 
(19), the final values of x and y at the termination of inte- 
eration, the number of changes of sign between successive 
addends in equation (16), the total number of calls made to 
SUBROUTINE VALUE during the process, and the clock time re- 


quired for numerical integration. 
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As an illustration consider the following example where 


ИЄ) ==> эз (21) 
s^ q^ 


A suitable form for SUBROUTINE VALUE is as follows: 


SUBROUTINE VALUE (X, Y, T, AL, U, V, MMMMM) 
e THIS IS TABLE ENTRY NUMBER 09 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL*8 (A,B,D-H,O-Z) 
IMPLICIT COMPLEX*16 (C) 
МММММ = MMMMM + 1 
ZERO=0.D+0 
CS=DCMPLX(X,Y) 


Note: The parameter a is 
CAL=DCMPLX (AL, ZERO) ( 
CT-DCMPLX(T, ZERO) d noted i by tie FORDRAN 


CDEN=CS**2+CAL**2 | 9) 

C=CS /CDEN 

CEXP=CDEXP(CS*CT) 

C=C*CEXP 

U=DREAL(C) 

V=DIMAG(C) 

RETURN 

END 
If the contour indicated by equation (18) were chosen and 
numerical inversion were performed upon the f(s) given by 
equation (21) for the case where A = 0.3, A = 0.1, a = 0.125, 
N = 3 and e = 1.D-11, the evaluation of F(27) would yield the 


following output 


АА F(2mT) |G(2T)| Хр Y. LOSC МММММ 


280 7.07107D-1 1.326350+0 -5.46 2.4 = 25 

where AA, LOSC, and MMMMM are the FORTRAN program names for 
the contour starting position, the number of sign changes 
between successive addends, and the total number of calls 
made to SUBROUTINE VALUE, respectively. Similar results 
for the almost 100 functions f(s) tested during the course 


of this research are listed in Appendix B. 
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IV. AUXILIARY SUBROUTINES 


Chapter III briefly mentioned the potential dangers in- 
herent in the utilization of complex arithmetic within SUB- 
ROUTINE VALUE when the function f(s) is multi-valued. 

If f(s) has branch points or winding points it is impor- 
tant to assure that one remains on the same branch of the 
function as one proceeds along the path of integration. 
Available FORTRAN operations do not assure this and it is 
necessary to employ specially programmed algorithms. The 
present discussion is limited to the two simple cases, rais- 
ing of a complex number to a non-integer power and taking 
the logarithm of a complex number. 

The essence of the matter lies in the FORTRAN operation 
DATAN2(y,x), or its equivalents, which must be relied upon 
to provide the argument of a complex number in polar form. 
This operation always provides a result in the range 
НО < т, SO that if Arg (2) actually passes through the 
values (2n+l)ı, for integer n, the FORTRAN produced value 
of Ars(z) experiences a jump of + 27. If one is determining 
the logarithm of Z or a non-integer power of z, shifting of 
Bench Will take place unless the continuity of Arg(z) is 
restored. 

SUBROUTINES CPOWER and CLOG provide for this continuity 
by monitoring changes in Arg(z) and adding or subtracting 


EX whenever it is appropriate to do so, to the value of 


22 





each provides 


As they are presenting written, 


to 
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for up to five separate calls, within SUBROUTINE VALUE, 
These subroutines are listed below, 


either SUBROUTINE CPOWER or SUBROUTINE CLOG. 


DATAN2(y,x). 
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K iLher algorithm is to be utilized by SUBROUTINE VALUE, 
the first requirement is to zero the applicable integer 
array, NCALL or MCALL, within the main program. Then, in 
SUBROUTINE VALUE, the calls to SUBROUTINE CPOWER are of 


ehe form 
CALL CPOWER(ZIN, ZOUT,P, J,NCALL) 


where ZIN is the COMPLEX*16 operand, ZOUT is the COMPLEX*16 
output, P is the REAL*8 power, and J is the INTEGER*4 index 
which tells whether this is the first, second,..., call to 
SUBROUTINE CPOWER within SUBROUTINE VALUE. The relationship 


obtained is 
р 
ZOUT = (ZIN) (22) 


SUBROUTINE CLOG is used Similarly. Calls made from SUBROUTINE 


VALUE are of the form 
CALL CLOG(ZIN, ZOUT, J, MCALL ) 


where 
ZOUT = Ln(ZIN) (2:32) 


J plays the same role as in SUBROUTINE CPOWER, There is no 


Р. It should be noted that the array MCALL replaces NCALL. 
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у. THE STEEPEST DESCENT CONTOUR 


The theoretical advantages of a steepest descent contour 
have been mentioned previously. Consider a contour in the 
complex plane along which v=constant. As shown in Figure 
2, let the € axis be tangent to the contour with the n axis 


perpendicular to the £ axis in the same sense that y 1S perpen- 


tentar to Xx. 


^t 
——————————————3» REAL AXIS 


Figure 2.  Steepest Descent Contour. Contour along 
which v is constant. Axis & is tangent at 
point P. Axis n is normal. 


Recall equation (11), 


g(s) - e9* £(s) = u(s) + iv(s) (11) 


By the Cauchy-Riemann relationships 


< 
Q2 
E 


ШО Чы _ dv _ du 

р St - $m on 
mene a vy = constant contour 

ЗУ O ou о (25) 

de ^ an 
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Such a v = с = constant contour divides the plane, locally 


Beast, into two regions: В. with v > ce and R, with v < с, 


1 2 
Consequently along such a contour =- is one signed, and, 
therefore, 37 is one signed. u + O along the v = c contour, 


in the direction of positive §. Therefore, there is no 
Ooscillation in u and, obviously also, there is none in v. 
meee u > 0 and v > O along a suitable integration contour, 
memeetection of c = 0 assures that this results. 

This second program to be discussed does partially follow 
a path of steepest descent. First it integrates vertically 
upward from a chosen point on the real axis (along a Bromwich 
contour) until intersection with a specified v = O contour 
occurs. Then it employs a tracking algorithm to obtain suc- 
cessive points on this contour along which it then integrates, 
moving to the left. 

SUBROUTINE VALUE has the same role as that described in 
Chapter III, namely, it calculates the u(s) and v(s) in 
accordance with equation (11) corresponding to a given x and 
Бот the initial part of the algorithm, the integration 
along the Bromwich contour, the value of x is constant, and 
y is increased in increments by Ay. Numerical integration 
in accordance with equation (16) is conducted just as it was 
done in the previous case. 

On the real axis, obviously, v = 0. With an initial 
incrementation in the vertical direction, v assumes a value 
Other than zero. The sign of v is recorded and the incrementa- 


tion process is continued until the sign of v changes. This 
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Benetitutes a contour crossing. The procedure continues 
until a predetermined number of such crossings are obtained, 
and then a simple linear interpolation is performed to locate 
Me specific ordinate value y for which v = O. Usually one 
might wish to use the first v = O contour above the real 
axis, but the second, third, or other such contour may be 
selected. 

The last (y; ,V;) values prior to the desired contour 
crossing and the first DETERE values after crossing are 
utilized in the linear interpolation process to predict YNEW, 
for which v = 0. SUBROUTINE VALUE provides the actual VNEW 
value corresponding to YNEW. Depending upon the sign of 
VNEW, the appropriate coordinate pair on i? Or En is 
replaced by (YNEW,VNEW) and the interpolation is repeated 
with this refinement. Termination of the interpolation occurs 
when YNEW results in a magnitude of VNEW which is less in 
magnitude than a prescribed tolerance limit. 

Once the values of x,y,u, and v have been established 
Lor the first point on the desired v = O contour in the pre- 
Vious process, the second, or tracking portion of the 
algorithm commences. 

SUBROUTINE XMARCH is called to provide a new coordinate 


Pair (x,y) in the negative sense of the real axis from the 


Known, or old, coordinate pair (х,у). This is accomplished 
by the following relationships 
х= х. — R cos(6) (26) 


О 


y = y, + R sin(9) 27) 








where (R,8) are the polar coordinates of a search pattern 
measured from the Known coordinates о) on the v = 0 
contour. The angle O is set at zero radians for the first 
call to SUBROUTINE XMARCH. SUBROUTINE VALUE is then called 

to provide the value of v at this point. The sign of v 
establishes upon which side of the v = O contour the new point 
lies. 

SUBROUTINE ROTATE is then called to increment the angle 
theta in an appropriate sense. А sweeping search is then 
conducted by SUBROUTINES ROTATE, XMARCH, AND VALUE until 
a Sign change is detected for v values between increments, 
їһи5,‚, indicating a crossing of the v = O contour. 

The last (9; ,v;) coordinate pair prior to a crossing and 
the first о coordinate pair after a crossing are 
utilized in a linear interpolation scheme, analogous to that 
of the first program segment, to refine the location of the 
new point on the v = O contour. 

Once the second point on the v = O contour is established, 
the tracking algorithm may be enhanced by predicting the 
third point on the contour to lie at the same angle 9 from 
the second point. As before, the sign of v at this new loca- 
tion permits a rotational search to be conducted in the 


appropriate sense. For small values of the radius vector 


R (R-.025 is built into the program, but may be changed by the 
user), this approximation has been found, in general, to be a 
reasonable one. 

The addends of equations (16) and 19) are evaluated at 


each step in the marching process within the main program 
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menee again, termination of the numerical integration 
occurs when each of N successive addend contributions is 
less, than a prescribed epsilon. 

Required inputs to the main program are the index speci- 
fying the desired v = O contour, the starting position on 
the real axis, the increment of the vertical march along the 
Bromwich contour, the termination criterion, any constants 
associated with the function f(s), and the value of t. These 
inputs are represented by the program variables NREG, AA, 
meer, EPS, (AL,BE,...), and T, respectively. 

The simple parameterized contour offers one decided ad- 
vantage over the steepest descent contour: each call made to 
SUBROUTINE VALUE is productive, in the former case, in the 
sense that it is utilized to compute a contribution of the 
addend of equation (16). The latter algorithm requires calls 
ESNSUSROUTINE VALUE to perform the linear interpolation pro- 
cedures, which are non-productive in terms of addend computa- 
tion. This represents increased computer time. The present 
investigation has shown that in each of the numerous cases 
which have been examined, the simple parabolic path indicated 


by equation (18) is more efficient and fully as accurate as 


a steepest descent contour. 








VI, SHANKS! ACCELERATOR 


The effectiveness of a family of non-linear sequence to 
Sequence transformations in accelerating the convergence of 
(some) slowly convergent sequences and in inducing conver- 
gence in (some) divergent sequences has been reported by 
Shanks [19]. If {A} (n=0,1,2,...) is a sequence of numbers, 
we form a Sequence of sequences [Ak ni which, for convenience, 
we write as {A(k,n)}. The integer k indicates the "order" 
of the sequence, with e - A(0,n); the integer n indicates 
the position of the term in the sequence. The rule for form- 


ing the sequence {A(k+l,n)} from the sequence {A(k,n)} is 
A(i,j) = 0 if j < 21 


| 2 
Me A(k,n)A(k,n-2) - A(k,n-1) (28) 
А(К,п)+А(К,п-2)-2А(К,п-1) 
Gt КЕ. ее 
To illustrate the power of this computational device con- 
Sider the application of equation (28) and its iterates to 


the first ten terms of the alternating series 


Ln(2) » 1 - C29» 


bol 
+ 

Cl 
| 

| 
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The results are 








п A(O,n) ACT) A) A(3,n) A(4,n) 
0 0.0 
.000сс00000 0.0 0-0 6: 
; b: 5696600000 0.0 0.0 0.0 9.9 
S 028233333333 0.7000009029 9.0 9.0 0.0 
4 С. 5832333333 2.690479190905 0.0 . my 2-2 
2 0: 8188666607 0.6924242424 0.693105753564 O. ae‘ 
2 (11545238095 0.6935891436 0.6931633407 0.693 1438893 0.0 
° Е. 99821355 98931808285 0: 5931473549 0.5931471951 
1d o 500006 9016936033417 0.6931451962 0.693 1471 119 056931471150 


where the tenth partial sum, A(0,10), is correct to only one 
Significant figure. However, iterate A(4,10) is correct to 
eight figures since Ln(2) = 0,6931471806. 

In this thesis research, Shanks' accelerator has been 
tested extensively in conjunction with the simple parameter- 
ized contours of integration discussed in Chapters II and 
BIT. 

Modification of the basic algorithm of Chapter III con- 
Sisted of the elimination from the main program of the ter- 
mination criterion for the numerical integration given by 
equation (16), and, instead, performing the integration for 
a finite number of terms, These addends of equation (16) 
were stored ina linear array and, subsequently, transformed, 
within the main program in accordance with equation (28). 

The accelerator did not enhance the rate of convergence 
of the numerical inversion process for any of the cases 
investigated. In fact, if the number of terms of А, which 
are transformed, becomes sufficiently large so that the 


Summation of the addend contributions of equation (16) 


approaches the correct result, the transformation will 








decrease the accuracy of the result as the denominator of 
equation (28) becomes small. 

Appendix A contains the program listing for the simple 
parameterized contour of integration in conjunction with 
Shanks' accelerator. 

Shanks [19] discusses transformations of higher order. 


These transformations may be worthy of investigation. 
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VII. LOCATION OF FUNCTION SINGULARITIES 


The criteria of Chapter II for the distortion of the 
Bromwich contour include the requirement that all of the 
Singularities of a function f(s) be to the left of the con- 
tour. Heretofore, the discussion has assumed that the loca- 
tion of these singularities was known. Indeed, this is not, 
in fact, always the case. In this chapter we investigate a 
method of determining the location of the singularities of 
f(s) when such information is either not known or is diffi- 
ore to Obtain analytically. 

To see how one can proceed without Knowing where the 


Singularities are, consider the simple case 


2 з ч. 
е 5 (30) 
S + o 
memechne particular values of a = 0.125 and t = Zt. 


We use a simple parabolic contour given by equation (13) 
with B=1 and with A, which we denote by the symbol хо in the 
discussion to emphasize its actual significance, being given 
various values. If we plot the result of the inversion 
versus the starting position Хо, the result is as shown in 
Eure 3. 

Por values of x less than zero, the result is 0.00. 

For Хо in the range from about 0.2 to 2.0, the result is 


07107. For Хо larger than about 2.0, the result is near 
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0.7 but is sensitive to the precise value of хо. Рог e in 
ШОШ се from 0.0 to 0.2, the result exhibits a "spike" 


Madea very significant change from 0.0 to 0.707107. 


Жо 
ze 
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Xo 
Eure 3. F(t) vs contour starting position for 
cos(at). Value of inverse of s 
fortc =F2 4S a function 2 2 
of starting abscissa x.. Anio) 


О 


We can explain this behavior as follows. For Xo less 
than zero, the singularities are on the wrong side of the 
contour and we simply get the wrong answer. For > in the 
neighborhood of 0.0 to 0.2, the contour passes so close to 
the singularity that the trapezoidal integration scheme is 
not accurate and experiences difficulty in getting an 
accurate value either for the right answer, 0.707107 or for 
the wrong answer 0.0. For Я Pnetherranse from 0,2 to 2,0 


the singularities are on the correct side of the contour and 








we get the correct answer. For Xo MAR AO ere are 
numerical difficulties, probably due to oscillation or to 

the effects of the exp(st) factor for s having large real 
part. 

We have done the same thing with many other transforms 
and have found the same behavior. Our experience leads to 
this suggestion for cases where the user does not know the 
location of the singularities. Perform the calculation for 
a number of different values of pe and plot the results as 
in Figure 3. The correct result is the ordinate of the 
rightmost plateau, before oscillations or other numerical 
difficulties have had a chance to introduce inaccuracies. 
In the case illustrated in Figure 3, there were only two 
plateaus, one corresponding to all the singularities lying 
on the wrong side of the contour and one corresponding to 
all the singularities lying on the correct side of the con- 
tour. However, with more singularities, one might expect 


to have several plateaus, only the rightmost of which repre- 


sents the correct result. 








VIII. THE HYPERBOLIC TANGENT 


In Chapter II it was stated that great advantage accrues 
from distorting the Bromwich contour as long as the process 
of distortion does not cause the contour to cross over a 
Singularity of f(s). In this chapter we investigate a typi- 
cal case in which it is impossible to satisfy this condition. 


The function 
E as 
f(s) = 5 tanh(->) (31) 


has an infinite sequence of Simple poles on the imaginary 
axis at points s = +(2n+1)— . By other methods it may be 
Shown that this function is the transform of the square wave 
function shown in Figure 4. In this chapter we investigate 
the extent to which our present methodology is capable of 


obtaining this square wave. 





Figure 4. Square wave function. 


The convergence of the numerical integration has been 


assured by bending the contour to the left in all other cases 








meperted in this thesis. Accordingly, we do so in the pres- 
ent instance. This implies that we must pass between two 
poles as the contour bends to the left and all of the in- 
finite number of poles above the point where the contour 
crosses the imaginary axis thus fall on the wrong side of 
the contour and are not included on the left. Thus, in 
effect, our procedure is finding the inverse Laplace trans- 
form of a substitute function which has only a finite number 
of poles coinciding in location with those of the given 
function f(s) and having the same residues. There is some 
reason to hope that the inverse of the substitute function 
will be essentially similar to that of the given function, 
at least for sufficiently small values of t. 

Indeed we do find this to be the case. SUBROUTINE CURVE 
of the algorithm of Chapter II was modified to formulate 
a simple parameterized curve, which utilized the Bromwich 
contour for numerical integration in its initial segment, 
starting from a position s=y on the positive real axis and 
continuing until the elevation of a finite number of poles 
is reached. The second segment of the simple parameterized 
curve which was used consists of the parabola of equation 
(18). This parabolic are allowed the integration contour 
to cross the imaginary axis between adjacent poles of the 
ШОО Оп f(s), Such a contour is indicated in Figure 5. 

The termination criterion which was used in this case 


is identical to that described in Chapter III. All other 


program input values remained unaltered from the form of their 
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Figure 5. Contour of integration for the hyperbolic 
tangent., Special contour for inversion of 
f(s) = 3 tanh(0s/2). It starts at x,, 
rises vertically, and then follows a para- 
bolic path between two poles. 
previous presentation with the exception that the integer 
JJJJJ, which selected the number of function poles to be 
encompassed by the contour was added. The program listing 
appears in Appendix A. 
The approximate solutions F(t) obtained by numerical 
inversion of the substitute function for the f(s) given 
by equation (31) are shown graphically below in Figures 6-9. 
Figure 6 is a reproduction of five periods of the square 
wave Of Figure 4 with a = 10.0 and the first 100 poles of 
f(s) above the real axis to the left of the contour of in- 
tegration. Figure 7 is a similar picture of the same F(t) 


with a = 5.0 and the first 150 poles of f(s) above the real 


axis encompassed by the contour, with only one period of the 


wave form. Figures 8 and 9 are analogous to Figures 6 and 7, 











Figure 6. Numerical approximation of five cycles of the 
square wave function witn 100 poles to the 
kerto Che contour. 
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Figure 7. Numerical approximation of one cycle of the 
square wave function with 150 poles to the 
left of the contour. 
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Iure 8. Numerical approximation of five cycles of 
the square wave function with 7 poles to 
the left of the contour. 


F (T) 





Figure 9. Numerical approximation of one cycle of the 
square wave function with 7 poles to the 
fete Or the contour. 








respectively; however, only the first seven poles of f(s) 
above the real axis were included on the correct side of the 
Esntour of integration. 

The reader should be aware that the increment of t used 
in the above series of graphs differed in each case. The 
inaccuracies shown in Figures 6-9 are partially due to the 
assumptions inherent within the employment of our algorithm, 
also, partially due to the fact that the plotter simply 
connects points sequentially. This can also be seen in the 
Боре ог the graphs at t = a, 2a, etc., which is a function 
of the particular At employed and not an inaccuracy of the 
Богі тит. 

Efforts to accomplish numerical inversion of the f(s) 
of equation (31) using the steepest descent contour were 
unsuccessful. The difficulty in using this contour appeared 
to be related to an adverse influence from the proximity of 
the function poles to the v = O contours near the imaginary 
axis. The algorithm was overtaxed in this vicinity and 
ШОШО not track the contour; consequently, further investi- 
gation was not attempted. 

In summary, we conclude that the effect of substituting 
a function which has only a finite number of poles which are 
coincident in location with the first n poles of a function 
f(s), which possesses an infinite number of poles spaced incre- 
mentally along the imaginary axis, produces an inverse 
approximation to the inverse of the given function f(s). 

The accuracy obtained with such a substitution is increased 
as more of the function singularities are encompassed by the 


distorted contour. 
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IX. FACTORS WHICH INFLUENCE ACCURACY 


Heretofore, our discussion has included the treatment of 
contour integration, utilizing a simple parameterized curve 
for that purpose, from two perspectives. The first of these 
methods, discussed in Chapter III, incorporated a termina- 
tion criterion which required a prescribed number, N, of 
successive integrand contributions each to have smaller mag- 
nitude than a specified epsilon. The second of these methods, 
discussed in Chapter VI, performed numerical integration for 
a finite number of points spaced along the contour with sub- 
sequent employment of Shanks' accelerator in an endeavor to 
enhance convergence. 

The first of these methods, we have found, was highly 
successful when employed with the proper point spacing, and 
a very tight termination criterion, namely, a small epsilon 
and large N. The second method has not, in our investiga- 
tions, been found to be successful, and need not be discussed 
further. The issue then becomes one of selecting the proper 
combination of step-size, epsilon, and N to provide an 
accurate and economical result. 

In order to make such an appropriate selection of these 
parameters, 1t is necessary to examine the factors which in- 


fluence accuracy. These factors include: the step-size along 


EsEEcontour, the termination criteria, and the oscillation. 








Consider a simple case, which we have chosen, along a 
path for which the oscillation is insignificant. The accu- 
racy in such a case is influenced, if the termination 
criteria are sufficiently stringent, only by the step-size 
along the contour of integration. Thus dealing with the 


transform 


pU e ш 


whose inverse is the cosine function, yields a plot of the 
logarithm of the percent error of F(t) versus the step-size, 
A, as shown in Figure 10. In this case, the contour of in- 
tegration was the simple parameterized curve of equation 
MES) with A = 0.3, B= 1.0, a = 0.125 and t = 27. 

The sequences of points in this plot are all obtained by 
maintaining N constant (N=3) and plotting 108, 0 percent error 
versus A with epsilon as a parameter. 

Four such parameteric plots appear in the figure below 
corresponding to the values of epsilon (EPS) as shown in 
the legend. 

Clearly the accuracy of the results may be seen to behave, 
in a general sense, in the manner shown in Figure 11, where 
EPS, < EPS, < EPS, 


When the termination criteria are tight, namely, when EPS 


< EPS . 


is small and N is, at least, greater than one, the percent 
error is a narrowly banded, roughly linear function of the 


integration step size, A, over a range from about A = 0.02 


to A = 0.2. This corresponds to EPS, in Figure ll. 
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Figure 10. Logarithm of percent error versus step size 
for the cosine function. 


If the termination criterion EPS is looser, as evidenced 
in other cases, shown schematically in Figure 11, the accu- 
racy of the result cannot be improved indefinitely by decreas- 
ing, A. Rather, there is a point, which differs in each case, 
beyond which the accuracy is decreased as the step size is 
decreased. The behavior is complicated by the fact that the 


effects of changes in step-size and changes in termination 


criteria are not themselves disjoint. Making the step-size 
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Figure ll. Schematic representation of accuracy factors. 
Behavior of the logarithm of percent error 
versus step size A as a function of the ter- 
mination criterion EPS, 

smaller results in meeting the termination criteria at an 
Parlier point along the path of integration and thus may 
actually reduce overall accuracy. 

One way of dealing with this matter, of course, is to 
revise the termination criteria so as to remove this inter- 
dependence. We simply have not experimented with using 
alternate forms for the termination algorithm. Instead, we 
Suggest the following viewpoint to the user who wishes to 
assure that he obtains results with a specified accuracy. 


He should first of all employ a very strict termination 


criterion and use a succession of rather small values of 








step size, and should also vary the starting point on the 
contour. In this way he should be able to establish an 
evaluation having even greater accuracy than he requires. 
Then, selecting what seems to be a good starting point and 
Keeping to it, he should increase the step size until the 
result is no more accurate than he requires. ше should 
observe how much saving in computer time (or in the number 
of calls to VALUE) this affords and should not increase step 
size beyond a reasonable point. Then he should loosen the 
termination criterion progressively, measuring the saving 
in computer time versus the possible loss of accuracy, 
stopping short of the point where he cannot rely on obtain- 
ing the required accuracy. 

In many cases where numerical integration is employed 
to obtain a result whose correct value is У, the numerically 


produced approximation y approximately satisfies a relation 


у = У + ал? (33) 


where a and m are unknown constants. The exponent m may 
frequently be established, once and for all, for the type 

of integration being used, and, for individual integrands, 
the correct value, Y, from two evaluations, y; obtained with 


А = А+ апа Yo obtained with A = A by using the extrapola- 


2? 


tıon formula 


Y = —— [à - А (34) 


оЎ1. 
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We have investigated this method of obtaining improved accu- 
racy for our inversions and have found that with the inter- 
dependence between termination criteria and step size, the 
numbers a and m were both rather large and difficult to 
establish - in other words, the extrapolation was not success- 
ful in producing improved values. Our concluSion was that 
it was more efficient from a computational point of view 
Simply to use a sufficiently small step Size and an appro- 
priately matched set of termination criteria so as to be 
able to obtain an accurate answer without interpolation. 
However, a suggestion for further study of this matter 


is made at the end of the next chapter. 
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Х. RESULTS, CONCLUSIONS, AND RECOMMENDATIONS 


This thesis has demonstrated that numerical integration 
along a suitable contour in the complex plane, thus implement- 
ing the complex inversion integral formula, is an effective 
way of obtaining numerical values for the inverse Laplace 
transformation of a given function. 

The success of the procedure, i.e., its efficiency, avoid- 
ance of inaccuracy due to oscillation, and its termination 
after including only a reasonable number of addends in the 
numerical sum, is related to the fact that the contours chosen 
bend to the left as they rise, thus taking advantage of the 
exp(st) factor in the integrand of the inversion formula. 

Although there are some theoretical advantages to what 
we have called the 'steepest descent" contour which follows 
a v=0 path, the algorithm which permits following such a path 
requires repeated evaluations of the integrand, many of which 
are not actually used in forming the summands for the numeri- 
cal summation. We have found by investigating almost 100 
cases for which known inverses are available in analytic form, 
that it is more efficient to use a simpler contour than the 
steepest descent path. By using what we have called a simple 
parameterized curve, we have devised a procedure which makes 


usSerof every evaluation of the integrand, We have found that 


the penalties of not following a steepest descent curve, 








namely an increase in oscillation of the integrand along the 
path, may be made quite tolerable by using any of a number 
of Simple curves. In particular, most of our successful work 
has been with a simple parabola. Our algorithm, however, 
includes te output of information which can serve to alert 
the user to the possibility of inaccuracy due to oscillation 
and thus suggest to him that he might select a different form 
for the path. 

We have even been successful in obtaining the inverse of 
a function which has an infinity of poles spaced along the 
imaginary axis. In doing so we violated the principle 
that all poles of the transform must be to the left of the 
contour. Nevertheless the results are quite satisfactory. 

Our success and accuracy have been so gratifying that 
we venture to suggest that our method may prove to be an 
efficient alternate method for the evalution of exotic func- 
tions for which other methods are slowly convergent or in- 
volve series the coefficients in which are difficult to 
obtain. For example, case 92 in Appendix B shows the suc- 
cessful evaluation of the rather uncommon Struve function. 

The original reason for attempting to employ numerical 
integration in the complex plane as a means of inverting a 
Laplace transformation was that alternate methods employed 
by Hiep and Zargary in conjugated heat transfer problems 
were Simply not accurate enough. They led to physically 


impossible results with some temperatures in the media below 


Sink temperature and others above source temperature. In 








this thesis we have not restudied these heat transfer prob- 
lems and recommend that this be done using our methods of 
inversion. There is some reason to believe that the diffi- 
culty encountered by Hiep and Zargary will not be encountered 
if our method is used, but we are not prepared to substan- 
Mate this claim at this time. 

Also it might be of use if the effects of varying the 
available parameters (shape of curve, spacing of points, 
Starting point of curves, termination criteria, etc.) were 
to be investigated further so that a user would be better 
prepared to deal with indications of inaccurate inversion 
or inordinate requirements for computer time. Our own 
numerical experiments were invariably so happily successful 
that we have not encountered need for such information. 

mer Tunctions f(s) for which it may be difficult to 
locate the singularities, we have shown that varying the 
Starting point of a Simple parameterized contour permits 
selecting an optimum contour in the sense that one can be 
eured that all singularities are to the left of the 
curve and also that the curve does not reach far enough to 
the right to impose numerical difficulties with large posi- 
tive exponentials. 

At the end of the preceding chapter we indicated that 
we did not find it profitable to employ extrapolation as a 
means of obtaining improved accuracy. However, this is 


probably worth looking at again, and our suggestions for 


doing so are as follows. First employ a termination criterion 








which is disjoint from the integration step size in the 

sense that the location, along the path of integration, of 
the point where termination occurs, is independent of step 
size. One way of doing this is to employ an epsilon in the 
termination criterion which is itself the result of multiply- 
ing a fixed, input epsilon by the step size. Then the 
extrapolation has a chance of being successful. So as to 
maintain optimum computational efficiency, one should use 
what is called Richardson extrapolation in which the evalua- 


tion points for the larger step size are themselves included 


among those for the smaller step size. 








APPENDIX A 


LISTING OF COMPUTER PROGRAMS AND SUBROUTINES 


The computer programs and principal subroutines that we 
have prepared and tested during our investigation are all 
listed within this appendix. Section A-1 contains the user's 
instructions glossary, flowchart, and program listing for the 
Simple parameterized curve. Section A-2 is a similar treat- 
ment for toe steepest descent contour. Section A-3 contains 
the user's instructions and program listing for the simple 
parameterized curve procedure adapted for use with Shanks' 
accelerator. Section A-4 contains the same information for 
adaptation of the simple parameterized curve for use with 
the nyperbolic tangent function described in Chapter VIII. 

It is hoped that sufficient detail has been provided to 
enable a user to adapt one of these programs to his purpose 
in an efficient manner. However, if additional insight is 
required, it may be necessary to refer to the chapter of 
this thesis in which the algorithm is developed. This is 
particularly true in the cases of Sections A-3 and A-4; the 
development in these sections has not been repeated where it 
is equivalent to that of Section А-1. 

Functions DREAL and DIMAG and SUBROUTINES CPOWER and 
CLOG, if they are required, are listed within Chapter III, 
as is a typical example of SUBROUTINE VALUE. Additionally, 
all functions f(s) which have been investigated during this 


research have their applicable SUBROUTINE VALUE listed within 


Appendix B. 








SECTION A-1: THE SIMPLE PARAMETERIZED CURVE 
USER INSTRUCTIONS 
This is the FORTRAN IV program for the numerical inver- 
Sion of the Laplace transformation of a function f(s), which 
performs contour integration along a distorted Bromwich con- 
tour in the form of a Simple parameterized curve. The fol- 
lowing instructions are provided to assist the user in 
familiarizing himself with the program so that he can adapt 
it to his particular requirements. 
1. SUBROUTINE CURVE is a user supplied group of instructions 
which formulates the numerical integration contour. The 
program listing presently contains a simple parabola, which 
may be utilized directly if desired. 
2. SUBROUTINE VALUE is a user supplied subroutine which 
calculates the real and imaginary components for g(s) = 
exp(st)*f(s) = u(s) + iv(s). For examples of preparation 
of this subroutine the user is referred to the test cases 
of Appendix B. 
3. The mandatory input variables to the main program are 
as follows: 
(а) AA - REAL*3 starting position (S=AA) on the real 
axis of the contour of integration. (Also denoted 
by A in equation (13) and Хо in Chapter VII.) 
(b) AL - REAL*8 value of any constant associated with 
the function f(s). (Note: If the function f(s) 
involves more than one constant, the calling state- 


ments to SUBROUTINE VALUE must be modified accord- 


iog ly: 








(с) DELP-REAL*8 increment of the p-parameter leading to 
spacing of points along the contour of integration. 

(d) Т - REAL*8 time for which F(t) is desired. 

4, Input parameters which may be altered by the user, at 
Bis discretion, are as follows: 

(a) EPS - REAL*8 tolerance parameter for the termina- 
tion of numerical integration. The program listing 
presently contains a value of 1.D-11. 

(b) NUMBER - INTEGER*4 number of successive addend con- 
tributions less than the specified value of EPS 
for which numerical integration is terminated. 

Ehe output variables which are printed as output by pro- 


gram are defined below: 


(a) AA starting position on the real axis 
(b) SUM the value of F(t) obtained by numerical 
M ersion ol f(s). (cf. equation (16)) 


(c) SUMA the absolute value of the summation of the 
the addend contributions of the numerical 
integracion AC. equation (14)) 

ed) X the final value of x on the contour at which 
integration was terminated 

(e) Y e inal value of y on the contour at which 
integration was terminated 

(f) LOSC the number of sign changes between successive 
addends encountered during numerical integra- 
tion 


(g) MMMMM the total number of calls made to SUBROUTINE 


VALUE during the numerical integration process. 








6. The output of these variables is as follows: 
AA, SUM, SUMA, X, Y, LOSC, MMMMM 
7. The time (seconds) required to perform the numerical inte- 
gration is printed on a line immediately preceding the items 
listed above. 
8. As previously described, the program will perform numeri- 
cal integration along one Simple parameterized curve. If the 
user desires, the program may be modified to "sweep" the con- 
Eur of integration over a range of starting values, AA, (cf. 
Chapter VII). This may be quickly accomplished in the follow- 
ing manner. 
(a) Replace the value of the integer KSTART from 1 to 
the number of contour integrations desired. 
(b) Replace the REAL*8 value of AA with the following 
statement 


AA = XSTART+AK*XINCRM 


where 
ASTART = the smallest starting value of tne 
contour of integration the user 
desires 
XINCRM = the Ax between starting positions 


of successive contours of integration 
must also be supplied,. 


Thus, the program will perform numerical inversion along 
ASRI contours of integration which differ only in their 
Starting positions along the x axis. The user may then 
observe increases or jumps in the value of the output variable 
SUM between successive integrations, enabling him to ascertain 
the locations of function poles. The utility of this proce- 


dure is discussed in depth in Chapter VII. 
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9, If SUBROUTINES CPOWER and/or CLOG are called by SUBROU- 
TINE VALUE in conjunction with use by this program, NCALL 
and/or MCALL must be dimensioned and initialized to zero 
within the main program and added to the calling statements 
to SUBROUTINE VALUE. Additionally, any powers, other than 


integer values, must be entered within the main program and 


passed within the calling statements to SUBROUTINE VALUE. 








VARIABLE 


AA 


ADD 


ADDA 


ADDDLD 


AL 


DELP 


EPS 


KOUNT 


KSTART 


LOSC 


MMMMM 


NN 


NUMBER 


PROD 





GLOSSARY OF VARIABLE NAMES 
DESCRIPTION 


Starting position on the real axis 
of the simple parameterized curve 


computed value of the addend of the 
numerical integration 


absolute value of the addend of the 
numerical integration 


saved value of the previous addend 
contribution 


constant associated with f(s) 


increment of point spacing along the 
Simple parameterized curve 


tolerance parameter for termination 
of the numerical integration process 


integer counter utilized to terminate 
the numerical integration process 


integer counter for the number of times 
the numerical integration is desired 

to be performed, using different start- 
ing positions 


integer counter utilized to record the 
number of sign changes between succes- 
sive addends 


integer counter to record the total 
number of calls made to SUBROUTINE VALUE 
during the numerical integration 


integer counter to record the total 
number of addend contributions 


integer input of the desired number of 
successive addend contributions less 
than a specified epsilon for termination 


parameter of the simple parameterized 
curve 


flag utilized to detect oscillation of 
the integrand 


LEGEND 





SUM total addend contribution for the 3 
numerical integration process 


SUMA absolute value of the total addend 3 
contribution for the numerical 
integration process 


T time for which F(t) is to be evaluated 2 


REST flag for termination of the numerical 1 
integration process 


U real part of the function g(s) 1 
UO previous value of U i 
V imaginary part of the function g(s) 1 
VO previous value of V 1 
Х real part of s 1 
XO previous value of X 1 
ү imaginary part of s 1 
YO previous value of Y 1 
LEGEND 


1. No action required by the user 
2. Mandatory user input 


3. Appears as program output 


4. May be altered by the user at his discretion 








START 


ESTABLISH THE 
SYSTEM LIMITS 
AND PARAMETERS 
INITIALIZE 
SET THE TIMER 


CALL 
CURVE AND VALUE 
OBTAIN 


INITIAL 
U AND V 





DO 400 
INCREMENT THE 
STARTING POINT 

COMMENCE 
NUMERICAL 
INTEGRATION 


99 


















COMPUTE THE 
INTEGRAND 
AND ITS 
ABSOLUTE 
VALUE 





SUM THE 
CONTRIBUTIONS 
OF THE 
INTEGRAND 
AND ITS 
ABSOLUTE VALUE 





60 


ESTABLISH 
MESS TEN ОЕ 
THE INTEGRAND 












COMPUTE THE 
PRODUCT OF THE 
OLD AND NEW 


INTEGRANDS 


LOSE 
TE THE COUNTER 
T.O.D+O FOR THE 
PRODUCT NUMBER OF 
SIGN CHANGES 
IS INCREMENTED 
KOUNT 
THE COUNTER 
LT.EPS | FOR INTEGRAND 
MAGNITUDE 
IS INCREMENTED 


LT.NUMBER 









GO TO 200 


KOUNT 








STOP THE 
TIMER 


WRITE 
THE 
OUTPUT 





400 
CONTINUE 


STOP 
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SECTION А-2: THE STEEPEST DESCENT CONTOUR 
USER INSTRUCTIONS 

This is the FORTRAN IV program for the numerical inver- 
sion of the Laplace transformation of a function f(s), which 
performs contour integration along the steepest descent path 
in the complex plane. The following instructions are pro- 
vided to assist the user in familiarizing himself with the 
program so that he can adapt it to his particular require- 
ments. 
1. SUBROUTINE XMARCH is a subroutine which tracks a speci- 
fied v=0 contour (steepest descent) in the complex plane, 
EC this contour has been intersected for the first time. 
As shown in the program listing, the following two equations 


X = XHOLD - 1.D-1*DCOS(THETA) 


Е YHOLD + 1.D—1*DSIN( THETA) 

march the numerical integrations in the favorable direction 
along the contour. The equations are parametric in THETA, 
which is the search angle provided by SUBROUTINE ROTATE. 

The 1.D-1 is a polar radius which may be altered to a larger 
or smaller increment by the user, or accepted as it appears. 
2. SUBROUTINE ROTATE increments the angular search angle 
THETA in order to locate a new point on the v-O contour. The 
angle is presently incremented by 5.D-2 radians with each 
call to the subroutine. This may be altered as required. 

3. SUBROUTINE VALUE is a user supplied subroutine which cal- 


culates the real and imaginary components for g(s) = exp(st) 


x Es) = u(s) + iv(s). For examples of preparation of this 
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Subroutine, the user is referred to Appendix B. 
4, The mandatory input variables to the main program are 
listed as follows: 

(a) AA - REAL*8 starting position (s=AA) on the real axis 

(b) AL - REAL*8 value of any constant associated with 
the function f(s). (Note: If the function f(s) 
involves more than one constant, the calling 
statements to SUBROUTINE VALUE must be modified 
Сеоне) 

(с) DELY - REAL*8 increment of the point spacing along 
the Bromwich contour. (Note: to alter the point 
Spacing along the v=O contour, the user is 
referred to 1 above.) 

(а) NREG - INTEGER*4 input value which specifies the 
first, second, third, etc. v=0 contour above 
the positive real axis, which is to be followed 
in the numerical integration process. 

(е) Т - REAL*8 time for which F(t) is to be evaluated. 

5. Input parameters which may be altered by the user, at 
his discretion, are listed below as follows: 

(a) EPS -— REAL*8 tolerance parameter for the termina- 
tion of numerical integration. The program 
listing presently contains a value of 1.D-ll. 

(b) NUMBER - INTEGER*4 number of successive addend 
contributions less than the specified value of 


EPS for which numerical integration is terminated. 
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аге 


The output variables which are received from the program 

defined below: 

(a) MMMMM - the number of calls made to SUBROUTINE VALUE 

(b) SUM -the total addend contribution for the numeri- 
cal integration 

(с) SUMA -the absolute value of the total integrand 
contribution for the numerical integration 

(d) EL - the computation time in seconds required for 
the process 

A sample program output is as follows 


694 MMMMM 


0.724724D+0 0.324972D+1 SUM, SUMA ко of vertical 
0.7071970+0 0.326734D+1 SUM,SUMA) 24 ¿na of 
0.71090D+1 EL computation 
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VARIABLE 


AA 


ADD 


ADDA 


AL 


CHECK 


DELY 


EPS 


EXAM 


KOUNT 


KSIGN 


LEVEL 


MMMMM 


NREG 


NUMBER 


GLOSSARY OF VARIABLE NAMES 
DESCRIPTION 


starting position on the real axis of 
the Bromwich contour segment of the 
integration 


computed value of the addend of numeri- 
cal integration 


absolute value of the addend of numeri- 
cal integration 


constant associated with f(s) 


flag to ascertain when a v=0 contour 
crossing has occurred 


increment of point spacing along the 
Bromwich contour segment of numerical 
integration 


tolerance parameter for termination of 
numerical integration 


flag to ascertain when a contour crossing 
has occurred 


integer counter used internally for test- 
ing contour regions 


integer counter utilized to terminate 
the numerical integration process 


integer counter utilized to record the 
number of v=0O contour levels crossed 


integer flag assigned to regions of 
positive and negative v 


integer counter for the number of calls 
made to SUBROUTINE VALUE during each 
segment of the numerical integration 
process 


the desired v=0 contour level above the 
positive real axis which is to be tracked 


integer input for the number of succes- 
Sive addend contributions less than the 
tolerance parameter required for termi- 
nation 
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LEGEND 


2 








SUM 


SUMA 


Т 


TERM 


TEST 


THETA 


THETAO 


THNEW 


THONE 


THTWO 


UHOLD 


VHOLD 


VNEW 


VOLD 


XHOLD 


YNEW 





СОЗИ ЧЕНО contribution for the 
numerical integration process 


absolute value of the total addend 
contribution over the contour of 
integration 

time for which F(t) is to be evaluated 


termination parameter for the absolute 
value of the addend 


Flag to ascertain when a contour 
crossing has occurred 


search angle used in the contour 
tracking algorithm 


previous value of THETA 
predicted angular location of a new 
point on the contour from a known 


previous point: 


coordinate utilized in linear inter- 
polation 


coordinate utilized in linear inter- 
polation 


пол вазе ое ==) 
previous value of U 
imasinary part of g(s) 
previous value of V 


actual value of the imaginary part of 
g(s) for the predicted v=0 location 


previous value of V 
real part of s 
previous value of X 


imaginary part of s 


predicted elevation of the desired v=0 
contour level 


K 








LEGEND 
No action required by the user 
Mandatory user input 


Appears as program output 


May be altered by the user at his discretion 
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START 


ESTABLISH THE 
PARAMETERS 
INITIALIZE 


Ser Tak 
TIMER 





CALL 
CURVE AND VALUE 

OBTAIN THE | 
INITIAL U AND V | 


200 
INCREMENT Y 
COMMENCE 
VERTICAL MARCH 


ql 








CALL 
CURVE AND VALUE 


COMPUTE THE 
VALUE OF THE 
INTEGRAND AND 
ITS ABSOLUTE 
VALUE 


SUM THE 
CONTRIBUTIONS 
OF THE 
INTEGRAND 
AND ITS 
ABSOLUTE VALUE 





SIGN 


EQ. 2 


ESTABLISH 


THE SIGN 
OF V 











СОМРПТЕ ТНЕ 
PRODUCT OF THE 
OLD AND NEW V 










K SIGN 
IF THE COUNTER 
LT.OD+O HOM Whee 
PRODUCT NUMBER OF 







CONTOUR LEVELS 
IS INCREMENTS 





300 


EEC NREG 


EONDUETZEINEAR 
INTERPOLATION 
LOCATE 


DESIRED POINT 
ON THE 
V=0 COUNTER 





3 








RETAIN THE 
Хх, О Гоу 
VALUES OF 


THE POINT 
ON THE COUNTER 


COMMENCE THE 
TRACKING ROUTINE 
TO LOCATE NEW 
POINTS ON THE 
COUNTER 


1100 
CALL 
X MARCH 
ROTATE 
AND 
VALUE 





ESTABLISH 
THE SIGN 


EQ. 2 


ОК у 
SET OLD V 











ЕХАМ 
COMPUTE THE 
PRODUCT OF 


OLD V AND V 





1100 


GT.O.D+0 LT.O.D+0 


1200 
1300 


CALL 
VALUE 
OBTAIN ACTUAL 
V NEW FOR THE 
INTERPOLATED 6 


COMPUTE 
HECK=V NEW * SIGN 








LT. EPS IF GUTTERS 
DABS(V) 








300 1300 





REPLACE THE 
COORDINATE 
PAIR ASSOCIATED 
WITH POSITIVE V 
REINTERPOLATE 






REPLACE THE 
COORDINATE 
PAIR ASSOCIATED 
WITH NEGATIVE V 
REINTERPOLATE 







1600 
COMPUTE THE 
INTEGRAND 
AND ITS 
ABSOLUTE VALUE 


SUM THE 


INTEGRAND 
AND ITS 
ABSOLUTE VALUE 












COMPUTE 
TERM=DABS( ADD) 







LT.EPS INCREMENT 
COUNTER 





1000 


ET. o EQ, 9 


WRITE OUTPUT 
STOP THE TIMER 
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LEGEND 


Mandatory user output. 


2. 


May be altered by the user at his discretion 


EST. 


may 


ОАТС into program, 


Polar radius R 
be changed by user 


9. 
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SECTION A-3.  SHANKS' ACCELERATOR 
USER INSTRUCTIONS 

The algorithm of this section is very similar to that 
of Section A-1. Only those features which differ have been 
discussed in this instruction. 
1. Contour integration is performed using a simple para- 
meterized curve for a finite number of points on this con- 
tour. This finite number of terms must be prescribed by 
the user. The integer variable M is used for this purpose. 
2. The dimension statements at the beginning of the program 
must provide storage space in linear arrays which are equal 
to or larger than the value assigned to M. (These quanti- 
ties will subsequently be used with Shanks! accelerator. 
Note particularly that the array is doubly indexed. ) 
3. The integer variable N specifies the number of iterates 
of Shanks' accelerator which are desired by the user. How- 
ever, the assigned value of N cannot be greater than 6 as 
the program is now written. 
4. A sample program output appears in Chapter VI of this 
thesis. Column 1 represents the natural result obtained 


after M evaluations of the addend, without application of 


the accelerator. 
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SECTION A-4. SPECIAL CONTOUR ALGORITHM FOR CASES AS 
TREATED IN CHAPTER EIGHT 
USER INSTRUCTIONS 

The integration contour used in this section is a vertical 
path to a specified value of y, followed by a simple parabola. 
The algorithm of this section is very similar to that of 
Section A-1. Only those items which differ are discussed in 
this section. 
1. Contour integration is performed in the initial segment 
of the program along the Bromwich contour. The integer vari- 
able JJJJJ specifies the number of poles of the function 
f(s) = š tanh(—) above the real axis that are to be enclosed 
to the left of the distorted Bromwich contour, which is used 
in the second segment of the program, 
2. The variable HALT computes the elevation above the real 
axis corresponding to the input value of JJJJJ for the 
specified hyperbolic tangent. 
3. The variable FLAG is internally assigned a value of zero 
or one, corresponding to the first or second segments of the 
contour of integration which the algorithm employs. No 
action is required by the user. 
4. The integer variable KSTART is utilized in the outer loop 
of the program to increment the value of t for which F(t) 
is computed. 
5. The linear arrays STORE and STASH are used to convert the 


values of F(t) and t to single precision, which is required 


for the graphics package associated with the IBM 360/67, 








6. SUBROUTINE VALUE contains the real arithmetic coding of 
me function g(s) for the hyperbolic tangent function of 
Chapter VIII of this thesis. Standard trigonometric identi- 
ties for the hyperbolic functions of a complex argument have 
been utilized. 

Note: This procedure is obviously keyed to the particular 
function f(s) = = tanh($). The user must modify it for 
other functions. One obvious modification which would 
generalize the procedure would simply be to specify rather 


than calculate the value of HALT, i.e., input HALT rather 


maon JJJJJ. 
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APPENDIX В 


LISTING OF TEST CASES AND EXPERIMENTAL RESULTS 


The Laplace transformation pairs f(s) and F(t) which we 
have used as test cases are listed within this appendix, 
Also included for each function f(s) is the listing of the 
SUBROUTINE VALUE and a typical result of numerical inversion 
at one specific value of t. Case 9 also shows a plot of 
F(t) versus t obtained numerically. 

The results obtained by numerical inversion have been 
compared where applicable, with tabulated values listed in 
(1). In all other cases the value of F(t) has been cal- 
culated on the IBM 360/67. 

The symbols which are used in the presentation of 


results are defined below: 


AL - constant associated with f(s) 

BE - constant associated with f(s) 

T - time for which F(t) is evaluated 

AA - real axis starting position of the integration contour 
X. - final value of x at termination 

Ye - final value of y at termination 


MMMMM - total number of calls made to SUBROUTINE VALUE 

LOSC - number of sign changes between successive addends 
SUM - value of F(t) (equation (16)) 

SUMA - value of G(t) (equation (19)) 

TIME - clock time (sec) required for numerical integration 


ANALYTICAL - analytical value of F(t) 


In all cases we used N = 5 and EPS = 1.D-11. 
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